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Abstract: We formulate a theory to show that the statistics of OCT signal 
amplitude and intensity are highly dependent on the sample reflectivity 
strength, motion, and noise power. Our theoretical and experimental results 
depict the lack of speckle amplitude and intensity contrasts to differentiate 
regions of motion from static areas. Two logarithmic intensity -based 
contrasts, logarithmic intensity variance (LOGIV) and differential 
logarithmic intensity variance (DLOGIV), are proposed for serving as 
surrogate markers for motion with enhanced sensitivity. Our findings 
demonstrate a good agreement between the theoretical and experimental 
results for logarithmic intensity-based contrasts. Logarithmic intensity- 
based motion and speckle-based contrast methods are validated and 
compared for in vivo human retinal vasculature visualization using high- 
speed swept-source optical coherence tomography (SS-OCT) at 1060 nm. 
The vasculature was identified as regions of motion by creating LOGIV and 
DLOGIV tomograms: multiple B -scans were collected of individual slices 
through the retina and the variance of logarithmic intensities and differences 
of logarithmic intensities were calculated. Both methods captured the small 
vessels and the meshwork of capillaries associated with the inner retina in 
en face images over 4 mm 2 in a normal subject. 
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1. Introduction 

The importance of retinal and choroidal vasculature visualization is paramount in diagnosing 
various eye diseases, such as diabetic retinopathy [1], age-related macular degeneration 
(AMD) [2], glaucoma [3], and anterior ischemic optic neuropathy (AION) [4]. 

Color fundus photography (CF) and fluorescein angiography (FA) have served as valuable 
imaging methods for retinal vasculature network visualization [5]. Deeper choroidal vessel 
imaging has also been realized by Indocyanine green angiography (ICG A) [5]. However, 
fluorescent dye injection can cause undesired side-effects [6,7], and the 2-D nature of these 
imaging techniques limits their applications for providing depth information and/or deep 
choroidal blood vessel visualization. 
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To meet the need for 3D retinal and choroidal vasculature assessment without the use of 
fluorescent dye injection, optical coherence tomography (OCT) has emerged as an attractive 
non-invasive depth-resolved imaging technology [8]. OCT realizes 3D visualization of retinal 
structure [9] and vasculature [10] by capturing interferograms of the light reflected back from 
the retina. Doppler OCT (D-OCT) [11] and Phase Contrast (PC)-OCT [12] are motion 
sensitive methods for vessel visualization, since these methods collect more than one depth 
scan of the same retinal loci and analyze Doppler phase shifts between depth scans at different 
depths. For example, a modified D-OCT method used autocorrelation and averaging 
techniques to capture retinal and choroidal vessels from complex OCT signals obtained in 
adjacent depth scans [13,14]. Although D-OCT captures the regions of high- velocity blood 
flow, such as in major vessels, it is unable to capture slow flow in retinal capillaries [11] or 
deep flows such as the choroidal circulation due to limited phase sensitivity and small time 
separation between two successive A-scans. Smart scanning protocols enhance sensitivity to 
the smaller signals expected from the microvasculature by increasing the time separation 
between two OCT depth scans and relying on the acquired phase [12,15-17] of OCT signals 
for contrast. However, these methods [15-17] such as PC-OCT [12] are highly sensitive to the 
phase instability of the system and environment by relying on the phase information. The 
required bulk motion removal [18] and timing-induced phase error correction algorithms [19] 
as well as extra optical module [19] add to the complexity of the phase sensitive swept source 
(SS)-OCT software and hardware. 

To avoid complications of phase instability in retinal microvasculature visualization, 
several elegant intensity contrast imaging methods have relied on the temporal speckle 
fluctuations caused by mobile particles, such as red blood cells. The first method captured 
microvasculature through Hilbert and Fourier analyses of the spectral domain (SD)-OCT 
signal and a moving reference arm [20,21]. However, the system complexity and high 
sensitivity to hyper-reflective static layers may limit its clinical application for distinguishing 
these areas from regions of motion. Another approach used speckle signals [22,23], which 
have been studied extensively in the literature [24-27]. Although this speckle variance 
imaging method was able to capture vessels no smaller than 25 um diameter size in the mouse 
dorsal skinfold [22,23], the multiplicative nature of the speckle made its variance highly 
sensitive to static regions similar to the first method [20,21], as shown in our results. 

Thus, there is a need for a simple OCT method which does not rely on the phase 
information and provides highly motion- sensitive contrast for distinguishing regions of 
motion from stationary areas. The latter is especially important for detecting leakage and 
abnormal vessels in patients with abnormal retinal and choroidal structure. 

Here, we purpose two motion- sensitive logarithmic intensity-based contrasts for 3D 
microvasculature visualization in an in vivo human retina. The proposed logarithm operation 
converts the multiplicative amplitude or intensity fluctuations (speckle) into the additive 
variations and recovers the motion contrasts by removing the speckle free signals (static 
regions) through statistical analysis. We develop two methods, called logarithmic intensity 
variance (LOGIV) and differential logarithmic intensity variance (DLOGIV), which provide 
these surrogate markers for motion by acquiring more than one depth scan of the same retinal 
voxel and analyzing logarithmic intensity fluctuations at the same voxel. To study and 
compare the proposed logarithmic contrasts with speckle amplitude and intensity contrasts, 
we formulate a theory for each. While theoretical and experimental results demonstrate the 
dependence of the speckle contrasts on the sample reflectivity strength and their dispensability 
for differentiating vessels from static regions, logarithmic contrasts show less sensitivity to the 
sample reflectivity strength. Application of LOGIV, DLOGIV, and speckle contrasts for 3D 
microvasculature imaging in the in vivo human retina is validated by employing a high-speed 
SS-OCT at 1060 nm. LOGIV and DLOGIV retinal en face views show the enhanced motion 
contrasts in comparison with speckle contrasts for capturing microvasculature that lies 
between hyper-reflective regions. Compared to the differential phase contrast method [12], 
these logarithmic intensity-based motion contrast methods are simpler, have similar 
performance, and do not require extra software and hardware [19]. 
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2. Theory 



The coherent summation of light waves with random path lengths causes destructive and 
constructive interferences called speckle [24]. The speckle pattern is stationary over time 
when the sample and optical imaging beam are static. However, backscattered light from 
moving particles, such as red blood cells, forms the speckle pattern that changes rapidly over 
time. Speckle has a multiplicative nature and its amplitude and intensity (OCT signal) follow 
Rayleigh and exponential distributions, respectively [25-27]. This highlights the potential 
application of speckle statistical analysis to capture blood vessels containing moving red 
blood cells [22,23]. However, the OCT signal sometimes contains specular components from 
strong (fixed) local scatterers, which drastically change the amplitude and intensity PDFs of 
speckle (OCT signal). In this scenario, amplitude and intensity statistics (such as mean and 
variance) depend not on only motion but also on hyper-reflective static regions due to the 
presence of these specular components, as shown in our formulated theory. Thus, amplitude 
and intensity-based (speckle) contrasts are not ideal for capturing and differentiating mobile 
scatterers from static regions. To remove the direct dependence of the speckle on the sample 
reflectivity, we test two different approaches: (1) mean scaling of the speckle contrasts 
(speckle contrast ratios) and (2) logarithmic intensity along with variance analysis. In order to 
study and compare these intensity-based contrasts, we develop a theory. Our theoretical 
results demonstrate high sensitivity of speckle contrast ratios to hyper-reflective static regions. 
However, the natural logarithm function and variance analysis separate the static signals 
(specular components) from speckle and cancel specular components (fixed terms), 
respectively. The developed theory shows that logarithmic intensity contrasts outperform both 
speckle contrasts and speckle contrast ratios in terms of motion sensitivity. 

Here, we first derive the generalized PDFs of the amplitude and intensity speckle in 
Fourier domain OCT. The Fourier transform of the received interferograms from a Fourier 
domain OCT system is intrinsically a complex stationary random process, and represents the 
amplitude (A) and phase {(/)) at a certain location and time as well as a zero mean complex 
additive white Gaussian wide-sense stationary random process (N(0,<7 n 2 )). 

U(z,t) = U Real (U) + jU lmag (z,t) = A(zj)e j ^ +N(z,t) (D 

where U Rea i and U Imag are the real and imaginary parts of the complex OCT signals, 
respectively. 

At a given depth (z) in the sample, the first term in Eq. (1) is proportional to a coherent 
sum of all backscattered light from N scatterers in a coherence volume with a size given by 
the probe beam area and coherence length of the light source. 

N 

U Re al (z, t) = A r (* 4 (x. , y. , z, 0 cos(^(x. , y . , z, 0) + N Re al (z, 0 = V x (z, 0 + N Re al (z, 0 ( 2 ) 

i=\ 

N 

U lmag (zj) = A r (0£ A s (x. , y. , z, t) sin(^(x. , y t , z, 0) + N lmag (z, 0 = V 2 (z, 0 + N lmag (z, 0 (3) 

i=l 

A n A s , N rea i(0,a nr 2 = a n 2 /2), and N imag (0,a n 2 = o 2 /2) are, respectively, the amplitude of light 
returned from the reference arm, the amplitude of backscattered light from a given scatterer in 
the sample arm located at (x h y h z), and zero mean real additive white Gaussian wide- sense 
stationary random processes. o n 2 and o n 2 are variances of the real and imaginary parts of a 
zero mean complex additive white Gaussian wide-sense stationary random process (N). In 
Eqs. (2) and (3), Vj(ju vl ,a vl 2 = o 2 /2) and V 2 (ju V 2, o v2 2 = o 2 /2) are the sum of many independent 
and identical random processes which approach two independent real Gaussian random 
processes with equal variances of a vl 2 and cr v2 2 , respectively. Thus, U(z,t) can be expressed as 
summation of a constant complex value (ju v i(z) + jjLi V 2(z)) and a zero mean complex Gaussian 
wide-sense stationary random process D(0, o 2 - o 2 + a 2 ). 
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U(z,t) = (ju vl (z) + jju v2 (z)) + D(z, t) (4) 
Linear intensity (I) is defined as 

I(z,t) = \U(z,tf =ULi(zJ) + Ul ag (z,t) (5) 

Thus, the intensity (I) and amplitude (Af = vf) of OCT signal have a non-central chi-square 
distribution with 2 degrees of freedom and a Ricean distribution [28], respectively. PDFs of / 
and Af are given by 

n(/) = 1 e^'^K Sl ^ 1 ) = —e °* I (2^L) ( 6a ) 

-(^ 2 +m 2 ) -(s 7 2 +m 2 ) 

where the non-centrality parameter of distributions and o 2 are 

s] = ju 2 vl + ju 2 2 (7a) 

^ = °nr 2 + ^ + ^vl 2 + ^v2 2 = ^C^/ + O ^ ) = g\ + C7* (7b) 

and /of J, the zeroth-order modified Bessel function of the first kind, is 

^^ir^+i) 2 

The non-centrality parameter shows the presence of specular components in OCT signal. 
Equations (6a) and (6b) approach exponential and Rayleigh distributions when the non- 
centrality parameter is zero (no specular component). While the thermal noise, shot noise, and 
relative intensity noise (RIN) of the laser source mainly determine o 2 , laser intensity 
fluctuation, transversal displacements between the scanning optical beam and sample, and 
mobile scatterers are contributing factors in cr v 2 

To find speckle amplitude variance (SAV), speckle intensity variance (SIV), speckle 
amplitude contrast ratio (SACR), speckle intensity contrast ratio (SICR), LOGIV, and 
DLOGIV as a function of OCT signal-to-noise ratio (SNR = s 2 /a 2 ), the mean and variance of 
the amplitude (Af), intensity (/), and natural logarithm of intensity (log(I)) are calculated. 

2.7 Speckle amplitude and intensity variances 

The mean and variance of the chi- square-distributed random variable with 2 degrees of 
freedom (/) and the Ricean random variable (Af) are [28] 

Ml = 2(a 2 nr + a\ ) + s 2 = a 2 + s 2 = cr 2 (1 + SNR) (9a) 
SIV = a 2 = 4(a 2 nr +a 2 vl ) 2 + 4(a 2 nr +a 2 Js 2 = a 4 +2a 2 s 2 =a 4 (l + 2SNR) (9b) 



1 a4n; <j4n 

2 l "-2«+o- v 2 1 ) / ~2~ °' 5 V r ' 2 

2 

SAV = a 2 M = 2(al+a 2 J + s 2 -^(a 2 nr+ a 2 vl )[L 0 ,( ^ )] 2 = a 2 {1 + SNR-^-[L 05 (-SNR)] 2 } 

(10b) 

where L 0 5 denotes a Laguerre polynomial. 
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Equations (9b) and (10b) clearly show that SIV (cr 2 7 ) and SAV (o 2 M ) depend not on only 
the noise power (<j n 2 ) and motion (cr v 2 ) but also on the sample reflectivity strength (s 2 ). Both 
SIV and SAV are also monotonically increasing functions of SNR. Therefore, speckle 
variance contrasts measure the noise power and motion of scatterers amplified with SNR 
value at the associated region. This highlights the inability of SIV and SAV to differentiate 
motion from static areas (cr v «0), especially hyper-reflective static regions with high SNR 
values. 

2.2 Speckle amplitude and intensity contrast ratios 

The speckle intensity and amplitude contrast ratios (SICR and SACR) are the ratios of the 
standard deviation to the mean of intensity and amplitude, respectively [29]. Using the derived 
Eqs. (9a) and (9b), and Eqs. (10a) and (10b), SICR and SACR can be expressed simply as 
following monotonically decreasing functions of SNR: 

SICR = ^J^SNR (na) 
jUj l + SNR 



SACR = ^= 4(l + SNR) -1 (Hb) 
Mm i^si-SNR)] 2 

Although these ratios cancel the direct dependence of speckle contrasts to the sample 
reflectivity, the derived Eqs. (11a) and (lib) demonstrate the SNR dependence of speckle 
contrast ratios and their mechanisms for capturing motion. Lower SNR is associated with 
higher speckle contrast ratios and motion. Monotonically decreasing nature of SICR and 
SACR makes these two contrasts less sensitive to hyper-reflective regions. However, their 
narrow dynamic ranges (1 and V47r _1 -1) and nonzero values of SICR and SACR at hyper- 
reflective static regions limit their application for capturing motion and differentiating it from 
structural signal, as described later in details. 

2.3 Logarithmic intensity and differential logarithmic intensity variances 

As described in appendix, the mean of the defined logarithm of a non-central chi-square 
random variable is 

iog(^)+r(o,^) = iog(^ 2 ) + JJ^ s t >o (12) 

V ' s I= 0 

\og(a 2 )-/ 

Therefore, the mean of the natural logarithm of the normalized intensity to the non-centrality 
parameter (s 2 ) is equal to T(0, s 2 /a 2 ). 

The variance of log(I), LOGIV, can also be given by 

n =\ n 

where Q is the normalized incomplete Gamma function [30] 

T(a) 

It follows that the variance of differences of logarithm intensities (two identical and 
independent random processes with equal variances) is 

Q(2n,SNR) 



DLOGIV = 2al g(I) =2 ^T J (15) 

n=\ n 
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LOGIV and DLOGIV approach 7i 2 /6 and 7i 2 /3, respectively, when SNR (s 2 /o 2 ) goes to zero 
(Eqs. (46) and (15)). 

To validate our theoretically derived Eqs. (12) and (13) for the mean and variance of the 
natural logarithm of intensity, one hundred thousand zero mean complex Gaussian random 
variables with a defined variance (a 2 ) were generated and added to the desired complex signal 
with different intensity values (ju vl 2 + ju v2 2 ), as expressed in Eq. (4). The red solid and blue 
dotted curves in Fig. 1(a) depict the derived Eq. (12) and the simulated result for the mean of 
the logarithm of the normalized intensity to the non-centrality parameter, r(0,(s/a) 2 ), as a 
function of different SNR values, respectively. The derived Eq. (12) clearly follows the 
simulated results. The mean of the logarithm of intensity approaches the logarithm of the non- 
centrality parameter for SNR>1. Figure 1(b) shows the derived Eq. (13) (blue dotted curve) 
and simulated result (red sold curve) for the variance of the logarithm of intensity. A close fit 
is found between the theoretical and simulated results in Fig. 1(b). The variance of the 
logarithm of intensity is 7r 2 /6 at a SNR equal to zero. When SNR is increased, the variance of 
the logarithm of intensity approaches zero. Thus, the motion sensing is based on the intensity 
fluctuation in logarithmic contrasts which are monotonically decreasing functions of SNR. A 
higher incomplete gamma function value detects lower SNR caused by the motion. 




SNR SNR 10log 10 (SNR)dB 



Fig. 1. The theoretical (blue dotted curves) and simulated (red solid curves) (a) mean of 
logarithm of normalized intensity to the non-centrality parameter and (b) variance of logarithm 
of intensity as a function of SNR. (c) The theoretical SACR, SICR, LOGIV, and DLOGIV as 
function of SNR. The close-up view shows the normalized SACR, SICR, LOGIV, and 
DLOGIV as function of SNR in Fig. 1(c). 

2.4 Comparison of speckle- and logarithmic intensity-based contrasts 

To capture motion caused by the intensity fluctuation of OCT signal, the ideal intensity-based 
motion contrasts need to be: (1) independent of the sample reflectivity strength (s 2 ) and (2) 
sensitive to motion over a wide SNR range. Theoretical results (Eqs. (9b) and (10b)) clearly 
show that speckle variances, SIV and SAV, depend on the sample reflectivity strength and are 
incapable of differentiating the motion signals from structural signals. 

Although mean scaling and logarithm operation along with variance analysis are able to 
remove the direct dependence of the speckle on the sample reflectivity (Eqs. (11a), (lib), 
(13), and (15)), speckle contrast ratios (SICR and SACR) and logarithmic-based contrasts 
(LOGIV and DLOGIV) vary with SNR (lla-b, 13, and 15). To compare SICR, SACR, 
LOGIV, and DLOGIV variations as a function of SNR, the derived Eqs. (11a), (lib), (13), 
and (15) were plotted in Fig. 1(c). The tail of normalized graph to the dynamic range 
describes the contrast sensitivity to hyper reflective-regions (Fig. 1(c), close-up view). While 
all contrasts decrease and approach zero with increasing SNR (Fig. 1(c)), logarithmic 
intensity-based contrasts demonstrate faster decline. Logarithmic intensity-based contrasts 
also show wider dynamic ranges than speckle contrast ratios (Fig. 1(c)). Comparisons of four 
normalized graph tails (Fig. 1(c), close-up view) demonstrate that LOGIV and DLOGIV are 
less sensitive to hyper-reflective (high SNR value) regions. 

Thus, less sensitivity to hyper-reflective regions make the logarithmic intensity-based 
motion contrasts more attractive than speckle contrasts for capturing mobile scatterers and 
differentiating motion from static regions in OCT applications such as ophthalmic OCT with 
SNR values < 35 dB. 



#159234 - $15.00 USD Received 1 Dec 2011; revised 31 Jan 2012; accepted 3 Feb 2012; published 10 Feb 2012 
(C) 2012 OSA 1 March 2012 / Vol. 3, No. 3 / BIOMEDICAL OPTICS EXPRESS 509 



3. Experimental methods 



3. 1 SS-OCT System 

To validate the proposed intensity-based methods for providing motion contrasts and compare 
them with each other and DPV method [12,17], we used a prototype 50 kHz phase sensitive 
SS-OCT system, incorporating a polygon-based 1060 nm (1015-1103) swept laser source 
[31], with -5.9 um axial resolution in tissue and 102 dB sensitivity (1.2 mW incident power). 
The SS-OCT system was comprised of the polygon-based swept-laser source [31], an 
interferometer, and a data acquisition (DAQ) unit (Fig. 2). The swept source output was 
coupled to the interferometer through an isolator where a 90/10 coupler was used to split light 
into a sample arm: reference arm. The sample arm light was split equally between the 
calibration arm and a slit lamp biomicroscope (Carl Zeiss Meditec) as shown in Fig. 2. A 
50/50 coupler combined and directed the reflected light from the sample to the one port of the 
interferometer output coupler. The reference arm light passed through a pair of collimators 
and was directed to the second port of the interferometer output coupler. The resulting 
interference fringes were detected on both output ports using a dual balanced photodetector. 
To generate a trigger signal at the beginning of the first interference fringe for data 
acquisition, a portion of the reference arm light was directed to a three port circulator and a 
fiber Bragg grating. A reflected optical pulse was detected using a photodetector and 
converted into an electronically tunable TTL signal as a trigger signal. The spectral signals 
were continuously digitized by triggering an AD conversion board. A D/A board was used to 
generate the driving signals of the two-axis galvanometers. A user interface and data 
acquisition was developed in Lab View to coordinate instrument control and enable user 
interaction. 



Fig. 2. Schematic of a polygon-based swept-laser source (gray box), interferometer (yellow 
box), and SS-OCT data acquisition unit used for imaging. SOA, semiconductor optical 
amplifier; diffraction grating (830 lines/mm); telescope (fi = 40 mm, f 2 = 35 mm); polygonal 
scanning mirrors (72 facets); COL, collimator; PC, polarization controller; ISO, isolator; CIRC, 
circulator; CUP, coupler; PD, photodetector; BD, dual balanced photodetector (InGaAs, 80 
MHz); FBG, fiber Bragg grating (0.1 nm). Data acquisition unit comprising of an AD 
conversion board (14-bit, GaGe CompuScope 14,200) for digitizing interference fringes and a 
D/A board (16 bit, National Instruments) for driving scanning mirrors (Cambridge 
Technology). 
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3.2 Scanning protocols 

The prototype SS-OCT instrument was used to image the eye of a healthy volunteer (3 7 -year- 
old Caucasian male) at the California Institute of Technology (Caltech). This study was 
approved by the Caltech committee for the protection of human subjects and adhered to the 
tenets of the Declaration of Helsinki. Written informed consent was obtained from subjects 
prior to OCT imaging. Total exposure time and incident exposure level were kept less than 5.5 
seconds and 1.2 mW in each imaging session, consistent with the safe exposure determined by 
American National Standards Institute (ANSI) and International Commission on Non-Ionizing 
Radiation Protection (ICNIRP). Patient interfaces were based on a Stratus OCT-3 system 
(Carl Zeiss Meditec) adapted with optics to support the 1060 nm wavelength range. A 60-D 
lens was used at the exit of the fundus camera with 13 mm working distance providing a beam 
diameter of 1.5 mm on the cornea (-15 um transverse resolution). 

Three scanning protocols were implemented. A ID protocol generated 1900 depth scans 
(A-scans) over the same transverse position in 0.038 seconds. A 2D protocol acquired four 
horizontal tomograms (B -scans) with 201 depth scans (A-scans) spanning the same transverse 
slice (2 mm) across the foveal centralis in 0.02 seconds. In the third protocol, a 3D OCT data 
set was collected by acquiring several neighboring B -scans over the parafovea. The system 
magnification, SS-OCT speed (50400 Hz), speed of the fast scan axis (200 Hz) with fly-back 
time (1 ms), and data acquisition time (4 seconds) gave an image size of 201 x 200 pixels over 
a 2 mm x 2 mm field of view (FOV); each B-scan was repeated four times. In the 3D 
scanning protocol, the fast scan axis was sagittal (superior-inferior) and the slow axis was 
horizontal (nasal- temporal). 

3.3 Image processing and motion contrast imaging 

The digitized signals were divided into individual spectral sweeps in the post-processing 
algorithm. Equal sample spacing in wave number (k) was achieved using a calibration trace at 
1.5 mm interferometer delay and numerical correction of the nonlinearly swept waveforms 
[32]. Image background subtraction and numeric compensation for second order dispersion 
[33] were performed. The complex SS-OCT data sets were upsampled by a factor of 4 and 
Fourier transformed. Axial motion correction was achieved on the obtained 2D and 3D SS- 
OCT data sets by cross correlating the consecutive horizontal tomograms. The estimated mean 
and variance of the linear intensity, and logarithm of intensity were calculated for all voxels 
through acquired depth scans. 

To perform motion contrast analysis and imaging, multiple A (or B)-scans were acquired 
over the same transverse position (or slice). Time separations were T A = 20 jus and T B = 5 ms 
between A- and B- scans for capturing the same position, respectively. Multiple linear 
intensity and phase measurements were recorded over the same transverse point separated in 
time. Six different intensity-based approaches were tested: SIV, SAV, SICR, SACR, LOGIV, 
and DLOGIV. In this paper, the number of measured samples (the linear intensity, logarithm 
of intensity, and phase) to estimate each of the different motion contrasts was 1900 in the ID 
protocol for each voxel through depth. In the 2D and 3D protocols, four samples were 
measured and used to estimate these motion contrasts. 

In the SICR, SACR, SIV, and SAV methods, the estimated amplitude and linear intensity 
means (/i), variances (a 2 ) as well as the ratios between their estimated standard deviations and 
means (cr//i) were calculated for the same transverse point acquired in successive A (or B)- 
scans (Figs. 3(b)-3(e)). LOGIV was realized by calculating the estimated variance of multiple 
logarithmic intensity measurements (LOG(I(z,T))) of the same transverse point acquired in 
successive A (or B)-scans separated in time (Fig. 3(f)). DLOGIV and DPV captured the 
differences between multiple logarithmic intensity (LOG(I(z,T))) and phase measurements 
((p(z,T)) of the same transverse points (separated in time) and calculated the estimated 
variance of these changes (Figs. 3(a) and 3(g)), respectively. To measure and remove timing- 
induced phase error [19] due to the random delay between the trigger signal and the 
subsequent A-to-D conversion (sample clock), a calibration signal was generated using a 
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stationary mirror in the calibration arm (Fig. 3) as described in [19]. The calibration signal 
was located at a depth of 2 mm in the OCT intensity image. The corrected phase differences 
between adjacent A (or B)-scans for the same transverse point at a given depth were 
calculated by subtracting the phase difference of the calibration signal, linearly scaled with the 
sample signal depth, from the measured phase differences (Fig. 3(a)) [19]. Phase unwrapping 
was performed on all measurements. A weighted mean algorithm [18] estimated and removed 
the bulk axial motion phase change error (Fig. 3(a)). 
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Fig. 3. Flowchart representing the data processing procedures for generating different motion 
contrasts of OCT: (a) DPV, (b) SICR, (c) SIV, (d) SACR, (e) SAV, (f) LOGIV, and (g) 
DLOGIV. 

The same described procedures were repeated for the adjacent transverse points in the 
same and neighboring B-scans to capture the retinal vasculature in 2D and 3D data sets. To 
remove SNR-limited intensity and phase change errors in 2D and 3D data sets for vasculature 
visualization, an average intensity threshold (10 dB above the mean value of the noise floor) 
was applied; all contrasts with average intensity values < mean (logi 0 (I no ise)) + 10 dB were set 
to zero in the corresponding images (Fig. 3). Additional processing was applied to the 3D 
contrast data sets using a moving average filter of ~5 um in the axial direction and a median 
filter of ~8 um radius in the transverse directions. 

To create the retinal en face views, the inner/outer photoreceptor segments (IS/OS) and 
vitreoretinal interface were detected using a segmentation algorithm [10]. Several depth 
integrated motion contrast en face images were generated by integrating the SICR, SACR, 
SIV, SAV, LOGIV, DLOGIV, and DPV between three different regions in the inner retina 
relative to IS/OS and vitreoretinal interface. 

4. Experimental results 

4.1 Experimental analysis of the linear intensity through retinal and choroidal layers 

In order to statistically analyze the linear intensity of the retina and choroid, a stationary 
optical beam was used to capture 1900 depth scans of a given transverse location through an 
in vivo human retina and choroid. The 2D foveal centralis tomogram (Fig. 4(a)) shows the 
transverse location and depths (red line) through which these A-scans were acquired in -38 
ms. Figure 4(b) presents the retinal and choroidal depth profile as a function of time (A- scan 
number) without alignment for the defined transverse position. This figure depicts the 
negligible axial bulk motion during collection of this data set. Dash-dotted blue and solid red 
lines show the estimated mean and standard deviation of the linear intensity through the retina 
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and choroid in Fig. 4(c), respectively. This figure clearly demonstrates SIV dependence on the 
layer reflectivity as described by Eq. (9b). The similar result can be found for the SAV. 

To study PDFs of the linear intensity and amplitude of OCT signal at different retinal and 
choroidal layers, the estimated non-centrality parameter (j/) and standard deviation of the 
measured linear intensity and amplitude (07 and a M ) were calculated throughout the captured 
depth profile. Figures 4(d)-4(o) show histograms of the measured linear intensity and 
amplitude at the level of six different retinal and choroidal layers, including the retina nerve 
fiber layer (RNFL), IS (inner segment), retina pigment epithelium (RPE), choriocapillaris 
(CC), Sattler's layer (SL), and Haller's layer (HL). The theoretically derived PDFs of the 
linear intensity (Eq. (6a)) are depicted as black dashed lines in Figs. 4(d)-4(i). These 
histograms follow chi-square distributions with 2 degrees of freedom (Eq. (6a)) for the 
estimated non-centrality parameter and standard deviation of the measured linear intensity. As 
shown in Eq. (6a), the PDF shape varies with the non-centrality parameter. For small non- 
centrality parameter values, the zeroth-order modified Bessel function of the first kind in 
Eq. (6a) is approximately equal to one. 
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Fig. 4. (a) Averaged intensity tomogram across the fovea centralis (5 mm). Red line depicts the 
targeted transverse location and depths by a stationary optical beam, (b) The retinal and 
choroidal structure was captured in 1900 repeated A-scans through depth (red line in Fig. 4(a)). 
(c) The estimated mean (dash-dotted blue line) and standard deviation (solid red line) of the 
linear intensity through the retinal and choroidal layers. Histograms and theoretically derived 
PDFs (black dashed lines) of the captured linear intensity at the level of (d) RNFL, (e) IS, (f) 
RPE, (g) CC, (h) SL, and (i) HL. Histograms and theoretically derived PDFs (black dashed 
lines) of the captured amplitude at the level of (j) RNFL, (k) IS, (l) RPE, (m) CC, (n) SL, and 
(o) HL. 
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In this scenario, the PDF of the linear intensity (Eq. (6a) approaches an exponential, as 
described previously [25-27]. Figures 4(d)-4(i) demonstrate a good agreement between 
experimental and theoretical results. Hyper- (Figs. 4(d)-4(f)) and hypo-reflective (Figs. 4(g)- 
4(i)) layers can be distinguished from each other using their linear intensity PDFs. Dashed 
lines in Figs. 4(j)-4(o) present the theoretically derived PDFs of the amplitude (Eq. (6b)). A 
close fit was found between histograms (Figs. 4(j)-4(o)) and a Ricean distribution (Eq. (6b)). 
The PDF of the amplitude (Eq. (6b)) matches with a Rayleigh distribution [25-27] when the 
non-centrality parameter approaches small values. 

4.2 Quantitative analysis of speckle contrast ratios and logarithmic intensity-based contrasts 
though the retinal and choroidal layers 

To investigate the quantitative variation of SICR, SACR, LOGIV, and DLOGIV through an in 
vivo human retina and choroid, a SS-OCT patient interface provided an optical beam for 
illuminating ~a 15 urn spot region on the retina of a healthy subject. The backscattered light 
was combined with the reference arm optical beam and detected in -38 ms. 1900 depth scans 
of a given transverse location through the retina and choroid were captured. SICR, SACR, 
LOGIV, and DLOGIV were calculated throughout 750 um depth by measuring the estimated 
mean (pi) and standard deviation (a) of the linear intensity and amplitude as well as the 
estimated variance of the logarithm of intensity and differences of the logarithm of intensities 
over the same data set as previously described in the motion contrast imaging section (3.3). 

Figures 5(a)-5(d) depict quantitative variations of SICR, SACR, LOGIV, and DLOGIV 
through the retina and choroid, respectively. The captured dynamic ranges of SICR, SACR, 
LOGIV, and DLOGIV are ~1, V(47r _1 -1), n 2 /6, and n 2 /3 as predicted by the theoretical and 
simulation results (Figs. 1(b) and 1(c)). Although LOGIV and DLOGIV values of hyper- 
reflective stationary layers such as RNFL, IS, and RPE are approximately zero (Figs. 5(c) and 
5(d)), the measured SICR and SACR values are nonzero for the corresponding regions 
(Figs. 5(a) and 5(b)). Nonzero values in Figs. 5(c) and 5(d) are associated with low SNR 
depths and/or regions of motion. These quantitative measurements (Figs. 5(a)-5(d)) are 
consistent with the theoretical results (Fig. 1(c)). LOGIV and DLOGIV demonstrate their 
potential application for capturing motion due to their wide dynamic range and low sensitivity 
to the sample reflectivity strength. Figure 5(c) shows that LOGIV values are ~7r 2 /6 in the 
choroid and sclera (>680 um). However, only DLOGIV exhibits different values over a wide 
range for different choroidal regions and sclera (Fig. 5(d)). This highlights DLOGIV 
capability for differentiating CC and SL from the outer choroid and capturing the choroid- 
sclera junction. 
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Fig. 5. Quantitative variations of (a) SICR, (b) SACR, (c) LOGIV, and (d) DLOGIV through 
the retina and choroid. 
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4.3 2D tomogram and en face view visualization of the retina using motion contrast imaging 
methods 

To study different motion contrast methods, four B-scans were acquired across the foveal 
centralis (2 mm). The averaged intensity of four obtained B-scans is depicted in Fig. 6(a). 2D 
SICR and SACR tomograms (Figs. 6(b) and 6(c)) show that these speckle contrast ratios 
capture not only regions of motion (between blue lines) in the inner choroid and small vessels 
(white arrows) in the inner retina but also highly reflective stationary regions in IS/OS, RPE 
(between red lines), and the inner retina. This finding is consistent with our theoretical result 
(Fig. 1(c)). While the SIV (Fig. 6(d)) is able to capture the inner retina vessels (white arrow), 
it highlights the static regions of IS/OS and RPE (between red lines) as motion. Motion in the 
inner choroid is barely detected in this tomogram. Figures 6(e) and 6(f) show the enhanced 
motion contrast in 2D LOGIV and DLOGIV tomograms. White static areas (between red 
lines) captured in 2D speckle tomograms (Figs. 6(b)-7(d)) are invisible in 2D LOGIV and 
DLOGIV tomograms (Figs. 6(e) and 6(f)). Regions of motion in the inner choroid (white band 
between blue lines) and the small vessels in the inner retina (white arrows) are detectable in 
these 2D tomograms (Figs. 6(e) and 6(f)). 

To compare the intensity-based contrasts with DPV contrast, 2D DPV tomograms are 
shown in Figs. 6(g) and 6(h) before and after compensation, respectively. Figures 6(g) 
demonstrate DPV is unable to capture motion without use of compensation algorithms [18,19] 
and an extra hardware module [19]. In addition, the calibration mirror image limits imaging 
depth [19]. Thus, the simplicity and motion sensitivity of LOGIV and DLOGIV may make 
these two contrast methods more attractive than other proposed phase- and intensity-based 
methods [12,15-17,20-23] for capturing motion and microvasculature. 




Fig. 6. Foveal (a) average intensity, (b) SICR, (c) SACR, (d) SAV, (e) LOGIV, (f) DLOGIV, 
(g) DPV before phase compensation, and (h) DPV after phase compensation tomograms (2 
mm). White regions correspond to regions with higher either motion or/and reflectivity. White 
arrows indicate the small vessel in Figs. 6(b)-6(h). IS/OS and RPE are signified between red 
lines (static regions). Regions of motion in the inner choroid are indicated between blue lines. 



Figures 7(a)-7(h) illustrate the inverted intensity, SICR, SIV, SACR, SAV, LOGIV, 
DLOGIV, and DPV en face views generated by integrating their values between the region 30 
um posterior to the vitreoretinal interface and the region 130 um anterior to IS/OS. Figure 7(a) 
shows that the meshwork of capillaries is barely visible in the intensity en face view. 
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Although small vessels and capillaries are seen in the SICR, SIV, SACR en face images 
(Figs. 7(b)-7(d)), the narrow dynamic range and high sensitivity to hyper-reflective static 
regions degrade retinal microvasculature enface visualization through contrast integration in 
the depth. Gray areas highlight the hyper-reflective stationary regions captured around the 
fovea avascular zone (FAZ) and between the interconnected microvasculature networks 
(Figs. 7(b)-7(d)). While SAV captures not only motion but also static regions (Fig. 6(d), 
between red lines), it is able to visualize the mesh work of capillaries through contrast 
summation in depth (Fig. 7(e)). Motion contrast enhancement is depicted in Figs. 7(f) and 7(g) 
using LOGIV and DLOGIV methods. Blood vessels in the ganglion cell layer and capillary 
meshwork of the inner plexiform layer are visualized in the LOGIV and DLOGIV en face 
views (Figs. 7(f) and 7(g)). FAZ is resolvable by considering the capillary network around it 
as shown in the LOGIV and DLOGIV images in Figs. 7(f) and 7(g). To compare retinal 
visualization using the proposed intensity-based motion contrast methods with the phase 
contrast method, the DPV en face image (Fig. 7(h)) is generated by summing DPVs over the 
same regions in the inner retina. Although LOGIV, DLOGIV, and DPV en face images 
(Figs. 7(f)-7(h)) achieve the similar contrast for foveal vasculature visualization, DPV is a 
complicated method due to its need for the compensation algorithms [18,19] and an extra 
optical module [19]. 




(f) (g) (h) 

Fig. 7. Parafoveal depth-integrated en face views over 4 mm 2 FOV acquired in 4 seconds. 
Inverted (a) averaged intensity, (b) SICR, (c) SIV, (d) SACR, (e) SAV, (f) LOGIV, (g) 
DLOGIV, and (h) DPV en face images of the inner retina. 
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To show the capillary meshwork of the inner retina through depth using logarithmic 
intensity-based motion contrast methods, the LOGIV and DLOGIV en face views are 
generated by integrating their values through different depths. Figures 8(a) and 8(b) show the 
capillary network of the inner retina between the regions 255 um and 216 um anterior to 
IS/OS in the inverted LOGIV, and DLOGIV en face views. The inverted DPV en face view 
(Fig. 8(c)) depicts the similar capillary meshwork of the inner retina in the same region. 
Similar retinal microvasculature network is also detected between the regions 216 um and 169 
jim anterior to IS/OS (Figs. 8(d)-8(f)) in the inverted LOGIV, DLOGIV, and DPV en face 
views. Figures 8(a)-8(f) clearly reveal depth-related variations of capillary meshwork 
morphology through the inner retina. 




Fig. 8. Parafoveal depth-integrated en face views over 4 mm 2 FOV acquired in 4 seconds. 
Inverted (a) LOGIV, (b) DLOGIV, and (c) DPV en face images of the retina between the 
regions 255 um and 216 um anterior to IS/OS. Inverted (d) LOGIV, (e) DLOGIV, and (f) DPV 
en face images of the retina between the regions 216 um and 169 um anterior to IS/OS. 



5. Discussion 

Our formulated theoretical description of the Fourier-domain OCT signal shows that the high 
sensitivity of amplitude and intensity speckle-based contrasts to static regions degrades 
microvasculature visualization. Inability of these contrasts to distinguish the motion signals 
from structural signals would also limit their application in detecting leakage and abnormal 
vessels in patients with irregular retinal and choroidal vessels and structure. To improve upon 
these methods, we described logarithmic intensity-based motion contrast methods and 
demonstrated their superiority for in vivo human parafoveal microvasculature visualization. 
Our results showed the potential application of depth-resolved LOGIV and DLOGIV methods 
for replacing invasive FA and ICGA for diagnosing eye diseases in the future. 
Microvasculature visualization in the retina with these logarithmic intensity-based motion 
contrasts also demonstrated their superior capability over previously reported phase [12,15- 
17] and intensity-based motion OCT [20-23] techniques without the need for compensation 
algorithms and extra optical modules. Further improvements to the logarithmic intensity- 
based motion contrast imaging methods might enhance vascular visualization, reduce the 
effects of inherent laser source intensity variation, and increase the scanning angle of the 
retina. 
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6. Appendix 

To calculate the mean and variance of the natural logarithm of the intensity (log (I)), we 
define a new chi-square random variable W (=I/o 2 ) with a non-centrality parameter of s w =s/o. 
Thus, log( I )can be written as 

log(7) = \og(a 2 W) = log(cr 2 ) + log(W) (16) 
where the PDF of W is given by 

p w (w) = e- (4+w) I 0 (2s w ^v) (17) 
The mean of the defined logarithm of a non-central chi-square random variable is 

m log(/) = log(a 2 ) + ]\og(w)e- (4+w) I 0 (2s w ^)dw ( 1 8) 

o 

To simplify Eq. (18), we use Eq. (8) and the n th derivative of the Gamma function given by 
[30] 



T (n) (k) = ^T(k) = J + V-V w [log(w)r dw (19) 



Therefore 



m log(/) =\og(cj 2 ) + e- sl f ^ r (1) (£ + l) (20) 

log(/) 5 tik\r(k+i) 

The first derivative of Gamma function can be expressed as [30] 

where yj 0 , y (=0.577), and H k (or H k (1) ) are the digamma function, Euler's constant, and 
harmonic number, respectively. 

By using the digamma function (21), the mean of log(I) in Eq. (20) can be rewritten as 

m^ ) =\o^)-y + e^± i ^^ (22) 

k=0 

The following identity Eq. (23) given by Gosper [34] simplifies Eq. (22) to Eq. (24) 

e -ff j (JtH± = r + mf) + log(f) (23) 

k=o 



=log(a 2 ) + log(4) + r(0,4) (24) 



where r is the incomplete gamma function. 
Since s w =s/a, the mean of log(I) is given by 



2 

io g (* 2 )+r(o,^) = iog(* 2 ) + [l^dt s t >o (25) 

* 1 s,=Q 

\og(a 2 )-/ 



Therefore, the mean of the natural logarithm of the normalized intensity to the non-centrality 
parameter (s 2 ) is equal to T(0, s 2 /a 2 ). 
The variance of log(I) is 
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<7 2 g(/) = E{ [\og(w)f } - [log(4 ) + r(0, 4 )] 2 (26) 

Using Eqs. (8), (17), (19) and the 2nd derivative of the Gamma function, the second moment 
of log (W) can be rewritten as 

E{ log 2 (w) } = e~ sl f (s2w)k r <2) (k + 1) (27) 

The second derivative of the Gamma function can be expressed as [30] 

l!!M = ^ + 1) + ^ (28) 

r(£+i) 0 dk 0 1 

where y/j(k+l) is the trigamma function and is related to the Riemann zeta function f and the 
generalized harmonic number of order k of 2 (H k (2) ) by [30] 

^(* + l) = £(2)-/^> = ~H™ =~t-\ (29) 

6 6 7=1 J 

Equations (21), (28), (29) and Taylor series expansion of exponential function simplify 
Eq. (27) to Eq. (30) 

£{io g 2 (w)} = ^ + r 2 - (2r)e~ 4 + ^ |;(4)W^ (3 o) 

6 k=o k\ k = 0 k\ 

By using the Gosper identity (23), the second moment of log(w) is rewritten as 

e{ log 2 ( w) } = ^ + y 2 - 2 r (r + iog(4 ) + r(o, 4 » + f(4 ) 0 1 ) 

6 



where F is defined by 



£=0 

To simplify Eq. (32), we define the first derivative of F as 



F(4) = g --X ( ^ )t[(g *| 2 " g " )] 02) 



^(g) = g - s y g -^"1 + tfn (33) 

where g= s w 2 and F(0)=0. 

Using the following recursion equations (34) and (35), the first derivative of F can be 
given by Eq. (36) 

H — = H (34) 

k+l k + \ k 

//£>-//< 2) =— ^ (35) 
The entire functions (E in ) [30] given in Eqs. (37), (38) simplify Eq. (36) to Eq. (39) for g>=0 
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g/i -I \ oo IT oo A: 

= ) (lz±A dt = = r+r(0, g ) + iog(g) = Z(-D w (f-) P7) 

£ . ( _ g) = _yll (38) 

^.Mz0^ w+ 2^ } } (39) 

g g dg g 

Therefore 

2 

F(4 ) = (y + io g (4 ) + r(o, 4 )) 2 + 2 f — (E in ( g ) + E in (- g ))d g m 

0 § 

Using Eqs. (26), (31), and (40), the variance of log(I) can be expressed as 

2 s w -g 

<(/) = — + 2S—(E in (g) + E in (-g))dg (41) 
6 o <? 

Now, we can simplify Eq. (41) using power series expansions of E in (Eqs. (37) and (38)) 



a 



2 _^ ^a + (-D*)K^4) 

log(/) 



2V v v ' " y ^^' (42) 



6 fc=1 kk\ 

where incomplete Gamma functions y(a,z) and are given by [30] 



^( fl ,z) = Jr fl -V^ ( 43 ) 



r(a,z) = j>-V'<fr 



(44) 



Using equalities given in Eq. (45) [30], the variance of /<?g(X) can rewritten as 

r(fl) = (fl - 1) != /(a, z) + T(a, Z ) (45) 
^(l + (-l) k )(T(k)-T(k,s 2 w )) 7T 2 *(! + (-!)*) *r(2n, 4) 



v ~ v _ ^ ^ ~ v ' ff// _ 2V I 4V h w ; (46) 



.2 

log(/) 



6 S ^! 6 ^ r S(2n)(2/i)! 

The following power series can be expressed as a Riemann zeta function [30] 



^ (l + (-l)*) = g(2) = ;r 2 
S ^ 2 2 12 



(47) 



and the variance of log(I) can be given by 



r(2«,4) 



^ =4 V °" ^f TQn,Sm (48) 

los<,) tT (2»)(2«)! ti (2«)(2»)! 



The variance of log(I), LOGIV, can also be given by 



WGN-- <(n --±®^ (49) 

n=l n 
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where Q is the normalized incomplete Gamma function [30] 

r(o) 

It follows that the variance of differences of logarithm intensities (two independent random 
processes with equal variances) is 



™ = 2< (I) =2£» (51) 
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